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Overview 
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•  Introduction 

•  Background 

•  Uncertainty  Analysis 

•  Systematic  Error 

•  Random  Error 

•  Conclusion 
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Goals 


w 


412TW 


■  Load  R,  JAGS  onto  your  laptop!  (Disk  set  up  for 
Windows) 

■  Learn  the  fundamentals  of  Bayesian  analyses 

■  Learn  how  to  run  Bayesian  analyses  from  within  R, 
using  JAGS,  and  interpret  the  results 

■  Learn  how  to  evaluate  “goodness  of  fit”  for  a  Bayes 
model 

■  Learn  predictive  posterior  distributions,  hierarchical 
modeling 


background 
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■  What  is  Bayesian  data  analysis?  Why  Bayes? 

■  Why  R  and  Bugs 

■  Bayesian  examples: 

■  binomial, 

■  normal  distribution 

■  reliability  applications 

■  Model  checking 

■  Bayes  estimate  and  prediction  of  lambda  in  HPP  reliability 
analysis 
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What  are  BUGS  and  R? 


V 


412TW 


•  Bugs  -  Bayesian  analysis  Using  Gibbs  Samplers 

■  BUGS  is  a  language  used  to  set  up  Bayesian  inference 

■  JAGS  (Just  Another  Gibbs  Sampler)  is  the  Bayes  software 
that  runs  within  R 

■  R  -  GNU  statistical  analysis  package 

■  Open  source  language  for  statistical  computing  and  graphics 

■  Well  vetted,  used  in  virtually  every  university  on  this  planet 

■  Bugs  from  within  R 

■  Offers  flexibility  in  data  manipulation  before  the  analysis  and 
display  of  inferences  after 

■  Avoids  tedious  issues  of  working  with  Bugs  directly 


A  brief  prehistory  of  Bayesian  data 

analysis 
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■  Reverend  Thomas  Bayes  (1763) 

■  Links  statistics  to  probability 

■  Laplace  (1800) 

■  Normal  distribution 

■  Many  applications,  including  census  [sampling  models] 

■  Gauss  (1800) 

■  Least  squares 

■  Applications  to  astronomy  [measurement  error  models] 

■  Keynes,  von  Neumann,  Savage  (1920’s-1950’s) 

■  Link  Bayesian  statistics  to  decision  theory 

■  Applied  statisticians  (1 950’s-1 970’s) 

■  Hierarchical  linear  models 

■  Applications  to  medical  trials,  conjugate  priors 

■  1990s  MCMC  techniques,  increased  computing  power 
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A  brief  history  of  Bayesian  data 
analysis,  BUGS,  and  R 


412TW 


■  “Empirical  Bayes”  (1950’s-1970’s) 

■  Estimate  prior  distributions  from  data 

■  Hierarchical  Bayes  (from  1970) 

■  Include  hyper  parameters  as  part  of  the  full  model 

■  Markov  chain  simulation  (from  1940’s  [physics]  and  1980’s  [statistics]) 

■  Computation  with  general  probability  models 

■  Iterative  algorithms  that  give  posterior  simulations  (not  point  estimates) 

■  R  code  (open  source)  for  statistical  applications  (1994) 

•  lme()  and  lmer()  functions  by  Doug  Bates  for  fitting  hierarchical  linear  and 
generalized  linear  models 

■  Bugs  (from  1994) 

■  Bayesian  inference  Using  Gibbs  Sampling.  Developed  explicitly  for 
Bayesian  statistics 
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What  is  Bayesian  data  analysis? 

Why  Bayes? 


■  Effective  and  flexible 

■  Combine  information  from  different  sources 

■  Examples  of  previous  uses  of  Bayes  from 
flight  test  include: 

■  Radar  systems  analysis 

■  Regression  testing  (“same  as  old”) 

■  Reliability  applications 

■  Multilevel  regression,  hierarchical  modeling-  test 
unit  parameters  are  not  all  the  same,  but  are 
drawn  from  “parent”  distribution 
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Structure  of  the  tutorial 


w 
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■  Computer  use 

■  Example  code  included  R,  and  JAGS 

■  “Follow  along”  computer  demonstrations 

■  Feel  free  to  Interrupt  with  questions 

■  Preliminaries  include 

■  How  to  set  up  a  BUGS  model  in  R 

■  Use  R  to  facilitate  posterior  distribution 
inference  and  diagnostics 
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Structure.,  continued 
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■  Understanding  how  BUGS  works  and  basic 

requirements  for  using  JAGS  with  BUGS 

■  Examples  -  Use  Bayesian  approach  to 

■  Estimate  the  parameter  of  a  binomial 
distribution 

■  Estimate  parameters  of  a  log-normal 
distribution 

■  Do  reliability  analysis-  examine  trend  in  an 
assumed  HPP 
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What  are  BUGS  and  R? 


V 
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■  BUGS  (Bayesian  Inference  Using  Gibbs  Sampling) 

■  Represent/Fit  Bayesian  statistical  models 

■  Is  a  “language”  designed  to  express  Bayesian  models 

■  R 

■  Open  source  language  for  statistical  computing  and  graphics 

■  BUGS  from  within  R 

■  Run  MCMC  based  Bayesian  analyses  from  within  R 

■  Offers  flexibility  in  data  manipulation  before  the  analysis  and 
display  of  inferences  after 

■  Avoids  tedious  issues  of  working  with  Bugs  directly 

■  [Open  R:  binomial] 
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Bayes,  Bugs,  and  R 


■ 
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■  Use  R  for  data  manipulations  and  various 
analysis  models 

■  Use  BUGS  within  R  to  fit  complex  Bayesian 
models 

■  User  R  to  summarize  results: 

■  Statistical  inference  from  a  posterior  distribution 

■  check  that  fitted  model  makes  sense  (validity  of 
the  BUGS)  result 

■  check  for  validity  of  model  implemented  in  BUGS 
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v  j  Fitting  a  Bayesian  model  in  R  and 
^  Bugs...  We’ll  cover 

■  What’s  required  for  a  BUGS  model 

■  Setting  up  data  and  initial  values  in  R 

■  Running  BUGS  and  checking  results 
(convergence,  model  adequacy) 

■  Displaying  the  posterior  distribution,  draw 
inferences 


412TW 
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V/  EVERY  R-script  using  JAGS  looks  like 

^ -  412TW 

1 .  Clear  the  workspace,  get  R2jags 
rm(list=ls()) 

require(R2jags)  #interface:  R  and  JAGS 

2.  Enter  the  “BUGS”  model  using  R-function 
cat()  ....  As  shown  on  next  slide 
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w  Class  Example:  Estimate  the  probability  of 
^  __juccess  of  a  rocket  launch  for  companies 

with  limited  launch/design  experience 
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•  Example  is  from  Hamada  et  al.,  Bayesian  Reliability, 
Springer,  2008 

•  Data:  1 1  companies  with  little  launch/design 
experience.  Objective  is  to  develop  a  statistical 
model  to  predict  launch  success  of  a  “new”  company 

•  Model  as  a  Bernoulli  process-  rocket  launch  was  a 
success  or  it  was  not 
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j  Example:  here  is  historical  data 


(1980-2000) 

412TW 

Vehicle 

Outcome 

Coded... 

Pegasus 

Success 

1 

Perch  eron 

Failure 

0 

AMROC 

Failure 

0 

Conestoga 

Failure 

0 

Ariane  1 

Success 

1 

India  SLV-3 

Failure 

0 

India  ASLV 

Failure 

0 

India  PSLV 

Failure 

0 

Shavit 

Success 

1 

Taepodong 

Failure 

0 

Brazil  VLS 

Failure 

0 
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Begin  with  Maximum  Likelihood 

_  Estimation  of  p 

■  Probability  of  success  is  p,  failure  is  (1-p) 

%ln.P)  =  (y)  Py(l  -  p)"-5’ 

■  Log-likelihood:  log[f(y|n,p)]  ay*  log(p)  - 
(n-y)  *  log(l-p) 

■  y  =  3,  n=1 1 ,  take  first  derivative  of  log- 
likelihood,  set  =0, 

-  0  =  d(  log(f(y|n,p))/d(p)  =  y/p  -  (n-y)/(1-p), 
solve  for  p 

-  p  =  3/1 1  =  0.272 
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Enter  the  BUGS  model 
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•  cat(‘  #  model  is  a  character  string 
model  { 

for(i  in  1:n)  { 

x[i]  ~  dbern(theta) 

} 

theta  ~  dbeta(1 ,1 )  #prior  on  theta 
}’ ,  #  end  of  BUGS  model 
file=“fileName.txt”)  #  end  of  cat() 
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R  and  Bugs  for  classical  inference 
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■  Estimate  the  parameter  of  a  binomial  distribution 
using  R  /  BUGS 

■  Displaying  the  results  in  R  or  rmarkdown 

■  Use  two  priors  for  the  analysis 

■  “vague”  prior-  uniform  across  (0,1) 

■  “informative”  prior-  p  around  0.3 
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Required  to  run  jags:  data 


fileNameData=list(x=c(rep(1 ,3),  rep(0,8)), 
n=1 1) 
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•  NOTE  data  must  be  a  list() 
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W  Required  to  run  jags:  inits 

fileNamelnits  =  function()  { 
list(theta  =  rbeta(1 ,1,1 ))} 
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•  NOTE  inits  must  be  a  function,  return  a  list 
(allows  for  multiple  MCMC  chains) 

•  Inits  can  be  a  NULL  function-  i.e.  let  JAGS 
pick  initial  values  of  parameters 
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Required  to  run  jags:  parameters  to  save 


fileNameParms  =  cftheta”) 

•  NOTE:  parameters  must  be  a  text 
“collection”  (vector)  of  variable  names 
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Summary,  so  far... 
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1 .  data  must  be  a  list 

2.  Inits  must  be  a  function 

3.  parameters  must  be  a  vector  of  text  name(s) 
of  the  variable(s)  we  want  to  examine,  use 
for  inference 
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Run  jags 


w 
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•  Call  to  jags:  (from  within  an  R  script) 
fileNameJags=jags( 

data=fileNameData, 
inits  =  fileNamelnits, 
parameters.to.save  =  fileNameParms, 
model.file=“fileName.txt”, 
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Run  jags  continued 


n.iter  =  2000, 
n.thin  =  1, 
n.burnin  =  1000, 
n. chains  =  4, 
DIC  =  TRUE) 


Notice  it’s  all  case  sensitive! 


Put  together.... 


412TW 


fileNameJags=jags( 

data  =  fileNameData, 

inits  =  fileNamelnits, 

parameters.to.save  =  fileNameParms, 

model.file  =  “fileName.txt”, 

n.iter  =  2000, 

n.thin  =  1 , 

n.burnin  =  500, 

n. chains  =  4, 

DIC  =  TRUE) 
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J  To  get  some  diagnostics,  and  a  plot: 

^ -  412TW 

■  fileNameJagsMC2  =  autojags(fileNameJags) 

■  attach.jags(fileNameJagsMC2) 

■  plot(density(theta)) 
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Now  you  try  it!  (exercisel.R) 


■ 
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■  Exercise:  set  up  and  run  the  binomial 
distribution-  estimate  theta,  get  a  posterior 
density  function  of  theta 

■  Use  3  successes  in  1 1  trials 

■  Uniform  prior  distribution  on  p 

■  Parameter  0,  plot  posterior  of  0 

■  Repeat  using  a  beta  distribution  for  prior  p, 
parameters  (alpha=2.24,  beta=2) 
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\,J  Overview  of  Bayesian  data  analysis 
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■  Decision  analysis  for  reliability 

■  Where  did  the  “prior  distribution”  come  from? 

■  Simulation-based  model  checking 
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Result  dependent  on  prior! 
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■  Different  priors  yielded  different  results! 

■  One  can  incorporate  prior  information  into 
analyses 

■  Prior  distributions  may  be  useful: 

■  Suppose  we  do  a  reliability  test  and  have  no 
failures  in  31 1  hours  -  what  can  we  say  about 
MTBF? 
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Decision  analysis  for  reliability 


w 
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■  Bayesian  inference 

■  Prior(0)  +  data  +  likelihood(data|0)  =  posterior^) 

■  Where  did  the  prior  distribution  come  from? 
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Prior  distribution 


V 
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-  Example  of  Bayesian  data  analysis 

■  Binomial 

-  Assume  a  beta  prior  for  p 

-  Incorporate  data  to  update  estimate  of  p,  MTBF 

-  On  the  disk-  binomial. R 

-  HPP  model 

-  Number  of  failures  proportional  to  interval  length 

-  Poisson  model 

■  On  the  disk-  poisson. R 

■  In  both  cases:  model  is  flexible- 

-  add  arbitrary  time  intervals,  new  data 
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More  on  Bayesian  inference 
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■  Allows  estimation  of  an  arbitrary  number  of 
parameters 

■  Summarizes  uncertainties  using  probability 

■  Combines  data  sources 

■  Model  is  testable 
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OK,  let’s,  estimate  p(successful 


mw 

■  MLE  has  excellent  “large  sample”  properties, 
but,  not  so  good  for  small  to  medium 
samples: 

■  large  sample  properties  of  MLE  do  not  pertain 
to  complicated  applications 

■  MLE  is  not  appropriate  for  hierarchical  models 

■  MLE  does  not  work  well  when  parameters  are 
close  to  boundary  of  the  parameter  space 

■  Deriving  analytic  expressions  is  difficult  in 
high-dimension  situations 

■  All  of  these  difficulties  are,  of  course, 
eliminated  in  Bayesian  estimation 
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Fundamentals  of  Bayesian  Inference 

412JW 

■  Frequentist  estimation  includes  a  confidence 
interval-  i.e.  an  interval  that  will  contain  the  true 
value  of  the  parameter  some  specified 
proportion  of  the  time  in  an  infinite  sequence  of 
repetitions  of  the  experiment 

■  Bayesian  estimation  combines  knowledge  of  the 
parameter  available  before  sample  data  are 
analyzed  with  information  gathered  during  an 
experiment 

■  Update  the  estimate  of  the  parameter 

■  Summarize  knowledge  of  the  parameter  using  a 
probability  density  function 
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Bayes  fundamentals  -  the  mechanism 


412TW 


pi&  I  >0 


fiy  I 

mi  J7) 


m(y)  =  f  fiy  \  0^p>i0^d0 


p(0|y)  is  the  posterior  density  of  0 

p(0)  is  the  prior  density  of  0 

m(y)  is  the  marginal  density  of  the  data,  and 

f(y|0)  is  the  sampling  density  of  the  data 
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And  the  parameter  of  interest,  8 
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E(0)  =  \df(d\y)d0 


Once  we  get  f(0|y)  we  can  estimate  any  density-related  parameter! 
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The  prior  distribution 
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In  the  launch  vehicle  example,  0  is  the 
parameter  of  interest,  the  probability  of 
success  of  a  launch 

Prior  information: 

-  Diffuse:  0  can  be  anywhere  in  the  interval  (0,1) 

-  Informative:  more  specific  information  about  0 
may  be  available-  past  history  indicates  that  0 
is  concentrated  near  0.4 

We  will  look  at  the  launch  problem  using  first 
the  “diffuse”  (aka  vague)  prior  and  then  the 
“informative”  prior 


38 


Priors 
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A  priori  we  take  all  values  in  the  interval  (0,1) 
to  be  equally  likely  for  0:  p(0)  =  1,  O<0<1 

OR  we  use  previous  experience  with  launch 
vehicles  to  assert  that  the  probably  of  a 
successful  launch  is  around  0.55,  and 
choose  for  the  prior  a  beta  distribution  with 
parameters  a=2.4,  and  3  =  2 

-  Mean  of  the  beta  distribution  is  a/(a+3  )  = 
0.545,  and 

-  Median  of  the  beta  distribution  is  (a-1)/(a+3-1) 
=  0.583 
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beta  density 


Beta  prior  for  p 
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Prior  distributions  for  p 
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Likelihood  function 
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•  Likelihood  =  bernouli,  so  result  is  either  1 ,  or 
0  with  probability  p,  repeated  1 1  times 

•  Can  use  a  single  likelihood,  binomial-  three 
successes  in  1 1  trials 
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Now  let’s  use  Bayes  rule  to  estimate  posterior 
distributions  of  the  parameter,  6 
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•  Bayes  Rule:  posterior  a  likelihood  *  prior 

•  Implement  this  in  the  “BUGS”  language 

•  Call  “jags”  to  develop  the  estimate  of  the 
posterior  distribution,  f(0  |y) 
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EVERY  R-script  to  use  JAGS  does 

the  following 
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■  Clear  the  workspace,  get  R2jags 
rm(list=ls()) 
require(R2jags) 

■  Enter  the  “BUGS”  model  using  R-function 
“cat()”  ....  As  shown  on  next  slide 
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Enter  the  BUGS  model 


•  cat(‘ 

model  { 

for(i  in  1  :n)  { 
x[i]  ~  dbern(theta) 

} 

theta  ~  dunif(1 ,1 )  #prior  on  theta 
}’ ,  #  end  of  BUGS  model 
file=“fileName.txt”)  #  end  of  cat() 
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Required  to  run  jags:  data 


w 
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■  fileNameData  =  list(x=c(1, 0,0,0, 1,0, 0,0,1, 

0,0),  n=1 1) 

■  NOTE  data  must  be  a  list() 
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Required  to  run  jags:  inits 
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fileNamelnits  =  function()  { 

list(theta  =  runif(1 ,0,1 )) 

} 


•  NOTE  “fileNamelnits”  must  be  a  function, 
return  a  list  (allows  for  multiple  MCMC 
chains) 
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Required  to  run  jags:  parameters  to  save 

•  fileNameParms  =  c(“theta”) 

•  NOTE:  fileNameParms  must  be  a  text 
collection  of  one  or  more  variable  names 


412TW 
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Quick  Check:  need  to  input  data,  inits, 
_  and  parameters  to  save 


412TW 


1 .  data  must  be  a  list 

2.  Inits  must  be  a  function 

3.  parameters  must  be  a  collection  of  text, 
naming  variables  we  want  to  examine 
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Run  jags 
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•  Call  to  jags: 

•  fileNameJags=jags( 

data=fileNameData, 

inits  =  fileNamelnits, 

para  meters,  to.  save  =  fileNameParms, 

model.file=“fileName.txt”, 
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Run  jags  continued 


n.iter  =  2000, 
n.thin  =  1 , 
n.burnin  =  1000, 
n. chains  =  4, 
DIC  =  TRUE) 


Notice  it’s  all  case  sensitive! 


Put  together.... 


w 


fileNameJags=jags( 
data  =  fileNameData, 
inits  =  fileNamelnits, 
parameters.to.save  =  fileNameParms, 
model.file  =  “fileName.txt”, 
n.iter  =  2000, 
n.thin  =  1, 
n.burnin  =  1000, 
n. chains  =  4, 

DIC  =  TRUE) 
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W  Get  some  diagnostics,  and  a  plot 

■  fileNameJagsMC2  = 
autojags(fileNameJags) 

■  attach.jags(fileNameJagsMC2) 

■  plot(density(theta)) 
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Now  you  try  it! 
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Exercise:  set  up  and  run  the  binomial 
distribution-  estimate  theta,  get  a  posterior 
density  function  of  theta 

■  Use  25  successes  in  289  trials 

■  Uniform  prior  distribution  on  p 

■  Parameter  0,  plot  posterior  of  0 
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Decision  analysis  for  reliability 


w 
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■  Bayesian  inference 

■  Prior(0)  +  data  +  likelihood(data|0)  =  posterior^) 

■  Where  did  the  prior  distribution  come  from? 
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Prior  distribution 
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■  Example  of  Bayesian  data  analysis 

■  HPP  model 

-  Number  of  failures  proportional  to  interval  length 

-  Poisson  model 

■  On  the  disk-  poisson. R 

■  Data  model 

-  Flexible:  arbitrary  time  intervals, 

■  Add  data  as  it  is  acquired 
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Types  of  prior  distributions 
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■  Two  traditional  extremes: 

■  Non-informative  priors 

■  Subjective  priors 

■  Problems  with  each  approach 

■  New  idea:  weakly  informative  priors 

■  Illustration  with  a  logistic  regression  example 
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Bayesian  inference-  reliability 
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■  Set  up  and  compute  model 

■  Use  data  at  hand;  update  as  more  data  becomes  available 

■  Inference  using  iterative  simulation  (Gibbs  sampler) 

■  Inference  for  quantities  of  interest 

■  Uncertainty  distribution  for  mean  time  between  failures 

■  Model  checking 

■  Do  inferences  make  sense? 

■  Compare  replicated  to  actual  data,  cross-validation 

■  Dispersed  model  validation  (“beta-testing”) 

■  Set  up  model  checking  in  the  HPP  program 
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Bayesian  inference  -  summary,  so  far 
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■  Set  up  and  compute  model 

■  Use  data  at  hand;  update  as  more  data  becomes  available 

■  Inference  using  iterative  simulation  (Gibbs  sampler) 

■  Inference  for  quantities  of  interest 

■  Uncertainty  dist  for  mean  time  between  failures 

■  Model  checking 

■  Do  inferences  make  sense? 

■  Compare  replicated  to  actual  data,  cross-validation 

■  Dispersed  model  validation  (“beta-testing”) 

•  Set  up  model  checking  in  the  HPP  program 
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Bayesian  inference 
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■  Allows  estimation  of  an  arbitrary  number  of 
parameters 

■  Summarizes  uncertainties  using  probability 

■  Combines  data  sources 

■  Model  is  testable  (falsifiable) 
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Model  checking 
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■  Basic  idea: 

■  Display  observed  data  (always  a  good  idea  anyway) 

■  Simulate  several  replicated  datasets  from  the 
estimated  model 

■  Display  the  replicated  datasets  and  compare  to  the 
observed  data 

■  Comparison  can  be  graphical  or  numerical 

■  Generalization  of  classical  methods: 

■  Hypothesis  testing 

■  Exploratory  data  analysis 

■  Crucial  “safety  valve”  in  Bayesian  data  analysis 
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Model  checking  and  model  comparison 
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■  Generalizing  classical  methods 

■  t  tests 

■  chi-squared  tests 

■  F-tests 

■  R2,  deviance,  AIC 

■  Use  estimation  rather  than  testing  where 
possible 

■  Posterior  predictive  checks  of  model  fit 

■  DIC  for  predictive  model  comparison 
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4  >  Model  checking:  posterior  predictive 


tests 


■  Test  statistic,  “T(y)” 

■  Replicated  datasets  y.rep(k),  k=1,...,n.sim 

■  Compare  T(y)  to  the  posterior  predictive 
distribution  of  T(y.rep(k)) 

■  Discrepancy  measure  T(y,theta(k)) 

■  Look  at  n.sim  values  of  the  difference,  T(y,thetak)  - 
T(y.repk,thetak) 

■  Compare  this  distribution  to  0 
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Model  comparison:  DIC  (deviance 


information  criterion) 


■  Generalization  of  “deviance”  in  classical  GLM 

■  DIC  is  estimated  error  of  out-of-sample 
predictions 

■  DIC  =  posterior  mean  of  deviance 

■  Compare  the  two  binomial  models: 

■  uniform  prior  (non-informative)  and 

■  beta(2.4,  2)  prior  (informative) 
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Understanding  the  Gibbs  sampler  and 

Metropolis  algorithm 
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■  Monitoring  convergence 

■  examples  of  good  and  bad  convergence 

■  n. chains:  at  least  2,  preferably  4 

■  Role  of  starting  points 

■  R-hat 

■  Less  than  1 .05  is  good 

■  Effective  sample  size 

■  At  least  1 00  is  good 
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Concluding  discussion 
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■  What  should  you  be  able  to  do? 

■  Set  up  hierarchical  models  in  Bugs 

■  Fit  them  and  display/understand  the  results  using 
R 

■  Compare  to  estimates  from  simpler  models 

■  Use  Bugs  flexibly  to  explore  models 

■  What  questions  do  you  have? 
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\J  Software  resources 
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-  Bugs 

■  User  manual  (in  Help  menu) 

-  Examples  volume  1  and  2  (in  Help  menu) 

■  Webpage  (  )  has  pointers 

to  many  more  examples  and  applications 

■  R 

-  ?command  for  quick  help  from  the  console 

■  Html  help  (in  Help  menu)  has  search  function 

-  Complete  manuals  (in  Help  menu) 

■  Webpage  (  )  has  pointers  to  more 

■  Appendix  C  from  “Bayesian  Data  Analysis,”  2nd  edition,  has 
more  examples  of  Bugs  and  R  programming  for  the  8-schools 
example 

-  “Data  Analysis  Using  Regression  and  Multilevel/Hierarchical 
Models”  has  lots  of  examples  of  Bugs  and  R. 
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